Comparison distance-based versus model-based approach
Packages
Data handling
Analysis
Custom functions
#--------------------------------------------------
# Plot ordinations from rda(), cca(), or cmdscale()
#--------------------------------------------------
plot_ord <- function(ord,comm, reverse_axis2 = FALSE){
site_scores <- scores(ord, display = "sites") %>%
as.data.frame() %>%
set_names("Axis1", "Axis2")
if(reverse_axis2){
site_scores$Axis2 <- -site_scores$Axis2
}
if(any(class(ord) %in% c("rda","cca"))){
sp_scores <- scores(ord, display = "sp") %>%
as.data.frame() %>%
set_names("Axis1", "Axis2")
}else{# PCoA
sp_scores <- wascores(site_scores, comm) %>%
as.data.frame()
}
ggplot() +
geom_point(data = site_scores, aes(x = Axis1, y = Axis2)) +
geom_path(data = site_scores, aes(x = Axis1, y = Axis2)) +
geom_text_repel(data = site_scores, aes(x = Axis1, y = Axis2, label = rownames(site_scores)), nudge_y = -0.1) +
geom_segment(data = sp_scores, aes(x = 0, y = 0, xend = Axis1, yend = Axis2), col = "grey") +
geom_text_repel(data = sp_scores, aes(x = Axis1, y = Axis2, label = rownames(sp_scores)), col = "grey") +
coord_fixed(ratio = 1) +
theme_minimal()
}
#--------------------------------------------------
# Plot ordinations from boral
#--------------------------------------------------
plot_boral <- function(boral_mod){
# Retrieve data to do the same plots as before
ord_data <- gg_lvsplot(boral_mod) %>%
layer_data(., 2) %>%
select(x, y, label) %>%
mutate(group = if_else(str_detect(label, "sp"), "species", "sites"))
ggplot() +
geom_point(data = ord_data %>% filter(group == "sites"), aes(x = x, y = y)) +
geom_path(data = ord_data %>% filter(group == "sites"), aes(x = x, y = y)) +
geom_text_repel(data = ord_data %>% filter(group == "sites"), aes(x = x, y = y, label = label), nudge_y = -0.1) +
geom_segment(data = ord_data %>% filter(group == "species"), aes(x = 0, y = 0, xend = x, yend = y), col = "grey") +
geom_text_repel(data = ord_data %>% filter(group == "species"), aes(x = x, y = y, label = label), col = "grey") +
coord_fixed(ratio = 1) +
theme_minimal() +
xlab("Latent variable 1") +
ylab("Latent variable 2")
}Legendre et Gallagher 2001
Create data of Fig. 3a
fig3a_data <- data.frame(
site = seq(1:19),
sp1 = c(7,4,2,1, rep(0,15)),
sp2 = c(1,2,4,7,8,7,4,2,1, rep(0,10)),
sp3 = c(rep(0,5),1,2,4,7,8,7,4,2,1,rep(0,5)),
sp4 = c(rep(0,10),1,2,4,7,8,7,4,2,1),
sp5 = c(rep(0,15), 1,2,4,7),
sp6 = c(0,1,1,rep(0,16)),
sp7 = c(rep(0,6),2,1,rep(0,11)),
sp8 = c(rep(0,11),3,1,rep(0,6)),
sp9 = c(rep(0,16), 4,1,0)
)
pos <- position_dodge(0.3)
fig3a_data %>%
gather(species, abundance, -site) %>%
ggplot(data= ., aes(x = site, y = abundance, fill = species, shape = species, group = species)) +
geom_line(aes(col = species), position = pos, alpha = 0.6) +
geom_point(data = . %>% filter(abundance != 0), alpha=0.6, size = 2.5, position = pos) +
scale_shape_manual(values=rep(c(21:25),4)) +
scale_x_continuous(breaks = seq(1:19)) +
theme_linedraw()Figure 1:
Reproduce Fig. 4
Using decostand :
- profile transformation: Y.tr = decostand(Y, “total”)
- chord transformation: Y.tr = decostand(Y, “norm”)
- log-chord transformation: Y.tr = decostand(log1p(Y), “norm”)
- Hellinger transformation: Y.tr = decostand(Y, “hellinger”)
- chi-square transformation: Y.tr = decostand(Y, “chi.sq”)
dist.ldc can also be used
comm <- fig3a_data
# Figure 4a : PCA on raw data
#-----------------------------
ord <- rda(comm)
fig_4a <- plot_ord(ord, comm) +
ggtitle("PCA on raw data")
# Figure 4b : CA on raw data
#---------------------------
ord <- cca(comm)
fig_4b <- plot_ord(ord, comm) +
ggtitle("CA on raw data")
# Figure 4c : PCA on chi2-transformed data
#-----------------------------------------
comm_tr <- decostand(comm,"chi.sq")
ord <- rda(comm_tr)
fig_4c <- plot_ord(ord, comm) +
ggtitle("PCA on chi2-transformed data")
# Figure 4d : PCA on chord-transformed data
#------------------------------------------
comm_tr <- decostand(comm,"norm")
ord <- rda(comm_tr)
fig_4d <- plot_ord(ord, comm) +
ggtitle("PCA on chord-transformed data")
# Figure 4e : PCA on species-profiles transformed data
#-----------------------------------------------------
comm_tr <- decostand(comm,"total")
ord <- rda(comm_tr)
fig_4e <- plot_ord(ord, comm) +
ggtitle("PCA on species-profiles transformed data")
# Figure 4f : PCA on Hellinger transformed data
#---------------------------------------------
comm_tr <- decostand(comm,"hellinger")
ord <- rda(comm_tr)
fig_4f <- plot_ord(ord, comm) +
ggtitle("PCA on Hellinger-transformed data")
# Figure 4g : PCoA with bray-curtis distance
#---------------------------------------------
comm_d <- vegdist(comm, "bray")
ord <- cmdscale(comm_d)
fig_4g <- plot_ord(ord, comm, reverse_axis2 = TRUE) +
ggtitle("PCoA with Bray-curtis")(fig_4a + fig_4b + fig_4c) /
(fig_4d + fig_4e + fig_4f) /
(fig_4g + plot_spacer() + plot_spacer()) +
plot_annotation(tag_levels = 'a', tag_prefix="(", tag_suffix = ")")Figure 2:
- Add variance explained to axes labels
- Do log-chord
Comparison of Fig. 4 with JSDM
Comparison with boral
Normal
# Normal - without row effects
#----------------------------
boral_normal <- boral(y = comm, family = "normal", lv.control = list(num.lv = 2), row.eff = "none")
boral_p1 <- plot_boral(boral_normal)+
scale_y_reverse() +
scale_x_reverse() +
ggtitle("Normal - without row effects")
# Normal - fixed row effects
#----------------------------
boral_normal_fixed <- boral(y = comm, family = "normal", lv.control = list(num.lv = 2), row.eff = "fixed")
boral_p2 <- plot_boral(boral_normal_fixed)+
scale_y_reverse() +
ggtitle("Normal - fixed row effects")
# Normal - random row effects
#----------------------------
boral_normal_rand <- boral(y = comm, family = "normal", lv.control = list(num.lv = 2), row.eff = "random")
boral_p3 <- plot_boral(boral_normal_rand)+
scale_y_reverse() +
scale_x_reverse() +
ggtitle("Normal - random row effects")
# Hellinger - Normal - without row effects
#----------------------------
comm_tr <- decostand(comm,"hellinger")
boral_normal_hell <- boral(y = comm_tr, family = "normal", lv.control = list(num.lv = 2), row.eff = "none")
boral_p4 <- plot_boral(boral_normal_hell)+
scale_y_reverse() +
scale_x_reverse() +
ggtitle("Hellinger-transformed - Normal - without row effects")
# Normal - fixed row effects
#----------------------------
boral_normal_hell_fixed <- boral(y = comm_tr, family = "normal", lv.control = list(num.lv = 2), row.eff = "fixed")
boral_p5 <- plot_boral(boral_normal_hell_fixed)+
scale_y_reverse() +
ggtitle("Hellinger-transformed - Normal - fixed row effects")
# normal - random row effects
#----------------------------
boral_normal_hell_rand <- boral(y = comm_tr, family = "normal", lv.control = list(num.lv = 2), row.eff = "random")
boral_p6 <- plot_boral(boral_normal_hell_rand)+
scale_y_reverse() +
scale_x_reverse() +
ggtitle("Hellinger-transformed - Normal - random row effects")(boral_p1 + boral_p4) / (boral_p2 + boral_p5) / (boral_p3 + boral_p6) +
plot_annotation(tag_levels = 'a', tag_prefix="(", tag_suffix = ")")Figure 3:
(boral_p1 + boral_p2 + boral_p3)+
plot_annotation(tag_levels = 'a', tag_prefix="(", tag_suffix = ")")Figure 4:
Poisson & Negative binomial
# Poisson - without row effects
#----------------------------
boral_poisson <- boral(y = comm, family = "poisson", lv.control = list(num.lv = 2), row.eff = "none")
boral_p1 <- plot_boral(boral_poisson)+
scale_y_reverse() +
scale_x_reverse() +
ggtitle("Poisson - without row effects")
# Poisson - fixed row effects
#----------------------------
boral_poisson_fixed <- boral(y = comm, family = "poisson", lv.control = list(num.lv = 2), row.eff = "fixed")
boral_p2 <- plot_boral(boral_poisson_fixed)+
scale_y_reverse() +
ggtitle("Poisson - fixed row effects")
# Poisson - random row effects
#----------------------------
boral_poisson_rand <- boral(y = comm, family = "poisson", lv.control = list(num.lv = 2), row.eff = "random")
boral_p3 <- plot_boral(boral_poisson_rand)+
scale_y_reverse() +
scale_x_reverse() +
ggtitle("Poisson - random row effects")
# Negative binomial - without row effects
#--------------------------------------
boral_negbin <- boral(y = comm, family = "negative.binomial", lv.control = list(num.lv = 2), row.eff = "none")
boral_p4 <- plot_boral(boral_negbin) +
scale_y_reverse() +
ggtitle("Negative binomial - without row effects")
# Negative binomial - fixed row effects
#--------------------------------------
boral_negbin_fixed <- boral(y = comm, family = "negative.binomial", lv.control = list(num.lv = 2), row.eff = "fixed")
boral_p5 <- plot_boral(boral_negbin_fixed) +
scale_y_reverse() +
ggtitle("Negative binomial - fixed row effects")
# Negative binomial - random row effects
#--------------------------------------
boral_negbin_rand<- boral(y = comm, family = "negative.binomial", lv.control = list(num.lv = 2), row.eff = "random")
boral_p6 <- plot_boral(boral_negbin_rand) +
scale_y_reverse() +
ggtitle("Negative binomial - random row effects")(boral_p1 + boral_p4) / (boral_p2 + boral_p5) / (boral_p3 + boral_p6) +
plot_annotation(tag_levels = 'a', tag_prefix="(", tag_suffix = ")")Figure 5:
Residuals of the models
- Normal - without row effects
## NULL
Figure 6:
Figure 7:
Figure 8:
Figure 9:
- Poisson - without row effects
## NULL
Figure 10:
Figure 11:
Figure 12:
Figure 13:
- Negative binomial - without row effects
## NULL
Figure 14:
Figure 15:
Figure 16:
Figure 17:
Comparison with Ecocopula ?
https://cran.r-project.org/web/packages/ecoCopula/vignettes/the_basics.html
Session info
## R version 3.5.0 (2018-04-23)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS 10.15.5
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] boral_1.8 coda_0.19-3 mvabund_4.0.1 vegan_2.5-6
## [5] lattice_0.20-35 permute_0.9-5 ggboral_0.1.7 patchwork_0.0.1
## [9] ggrepel_0.8.2 forcats_0.3.0 stringr_1.4.0 dplyr_1.0.2
## [13] purrr_0.3.2 readr_1.1.1 tidyr_1.0.0 tibble_2.1.3
## [17] ggplot2_3.3.0 tidyverse_1.2.1
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-137 lubridate_1.7.4 httr_1.4.1 tools_3.5.0
## [5] backports_1.1.4 R6_2.4.1 R2WinBUGS_2.1-21 mgcv_1.8-23
## [9] questionr_0.7.0 colorspace_1.4-1 withr_2.1.2 tidyselect_1.1.0
## [13] compiler_3.5.0 cli_2.0.2 rvest_0.3.4 xml2_1.2.2
## [17] labeling_0.3 bookdown_0.13 scales_1.0.0 mvtnorm_1.1-0
## [21] digest_0.6.25 rmarkdown_1.15 pkgconfig_2.0.3 htmltools_0.3.6
## [25] highr_0.7 rlang_0.4.8 readxl_1.1.0 rstudioapi_0.9.0
## [29] shiny_1.3.2 generics_0.0.2 jsonlite_1.6 magrittr_1.5
## [33] Matrix_1.2-14 Rcpp_1.0.5 munsell_0.5.0 fansi_0.4.1
## [37] abind_1.4-5 lifecycle_0.2.0 stringi_1.4.6 yaml_2.2.0
## [41] MASS_7.3-51.4 plyr_1.8.4 fishMod_0.29 grid_3.5.0
## [45] parallel_3.5.0 promises_1.0.1 crayon_1.3.4 miniUI_0.1.1.1
## [49] haven_2.1.1 hms_0.4.2 knitr_1.24 pillar_1.4.2
## [53] boot_1.3-20 corpcor_1.6.9 reshape2_1.4.3 codetools_0.2-15
## [57] R2jags_0.5-7 glue_1.4.2 evaluate_0.14 modelr_0.1.5
## [61] vctrs_0.3.5 rmdformats_0.3.5 httpuv_1.5.2 cellranger_1.1.0
## [65] gtable_0.3.0 assertthat_0.2.1 xfun_0.9 mime_0.7
## [69] tweedie_2.3.2 xtable_1.8-4 broom_0.5.2 later_0.8.0
## [73] rjags_4-10 cluster_2.0.7-1 statmod_1.4.34 ellipsis_0.3.0